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Abstract 

Effects of immune delay on symmetric dynamics are investigated within a model of antigenic 
variation in malaria. Using isotypic decomposition of the phase space, stability problem is 
reduced to the analysis of a cubic transcendental equation for the eigenvalues. This allows one 
to identify periodic solutions with different symmetries arising at a Hopf bifurcation. In the 
case of small immune delay, the boundary of the Hopf bifurcation is found in a closed form in 
terms of system parameters. For arbitrary values of the time delay, general expressions for the 
critical time delay are found, which indicate bifurcation to an odd or even periodic solution. 
Numerical simulations of the full system are performed to illustrate different types of dynamical 
behaviour. The results of this analysis are quite generic and can be used to study within-host 
dynamics of many infectious diseases. 

1 Introduction 

Among various strategies employed by pathogens to evade the host immune system, a prominent 
place is occupied by antigenic variation. Notable examples of pathogens relying on this strategy 
of immune escape include African Trypanosoma, Plasmodium falciparum, HIV, several members 
of Neisseria family, Haemophilus influenzae etc. [IB]. Despite the fact that some details of this 
mechanism are still unclear, the main features of this method of immune evasion are quite universal. 
The immune system of the host detects potential infection with a pathogen by identifying specific 
chemical determinants, such as proteins and carbohydrates, known as epitopes, on the surfaces 
of infected cells. This triggers the differentiation of precursor cells into effector cells, which are 
then able to eliminate the infection. Some pathogens have evolved to have a wide variety of surface 
markers (antigens), and by changing the antigens the present on the cell surface, these pathogens can 
for a long period of time remain unrecognized by the immune system, giving them an opportunity 
to be transmitted to other hosts. This process of sequentially presenting different antigens in order 
to avoid the host immune system is known as antigenic variation. 

There are several particular ways of implementing antigenic variation. In the case of Try- 
panosoma brucei, the organism that causes sleeping sickness, parasite covers itself with a dense 
homogeneous coat of variant surface glycoprotein (VSG). Genome of T. brucei has over 1000 genes 
that control the expression of VSG protein, and switching between them provides the mechanism 
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of antigenic variation [M]. What makes T. brucei unique is the fact that unlike other pathogens, 
whose antigenic variation is typically mediated by DNA rearrangements or transcriptional regula- 
tion, activation of VSGs requires recombination of VSG genes into an expression site (ES), which 
consists of a single vsg gene flanked by an upstream array of 70 base pair repeats and expression 
site associated genes (ESAGs). T. brucei expresses one VSG at any given time, and the active VSG 
can either be selected by activation of a previously silent ES (and there are up to 20 ES sites), or 
by recombination of a VSG sequence into the active ES. The precise mechanism of VSG switching 
has not been completely identified yet, but it has been suggested that the ordered appearance of 
different VSG variants is controlled by differential activation rates and density-dependent parasite 
differentiation [3U [51] . 

For the malaria agent P. falciparum, the main target of immune response is Plasmodium falci- 
parum erythrocyte membrane protein-1 (PfEMPl), which is expressed from a diverse family of var 
genes, and each parasite genome contains approximately 60 var genes encoding different PfEMPl 
variants [21 . The var genes are expressed sequentially in a mutually exclusive manner, and this 
switching between expression of different var gene leads to the presentation of different variant sur- 
face antigens (VSA) on the surface of infected erythrocyte, thus providing a mechanism of antigenic 
variation [10^141) . In all cases of antigenic variation, host immune system has to go through a large 
repertoire of antigenic variants, and this provides parasites with enough time to get transmitted to 
another host or cause a subsequent infection with a different antigenic variant in the same host. 

Despite individual differences in the molecular implementation of antigenic variation, such as, 
gene conversion, site-specific DNA inversions, hypermutation etc., there are several features com- 
mon to the dynamics of antigenic variation in all pathogens |16| . These include ordered and often 
sequential appearance of parasitemia peaks corresponding to different antigenic variants, as well 
as certain degree of cross-reactivity. Several mathematical models have been put forward that aim 
to explain various aspects of antigenic variation. Agur et al. [2j have studied a model of antigenic 
variation of African trypanosomes which suggests that sequential appearance of different antigenic 
variants can be explained by fitness differences between single- and double-expressors - antigenic 
variants that express one or two VSGs. However, this idea is not supported by the experimental ev- 
idence arising from normal in vivo growth and reduced immunogenicity of artificially created double 
expressors [40]. Frank [20] has suggested a model that highlights the importance of cross-reactivity 
between antigenic variants in facilitating optimal switching pattern that provides sequential dom- 
inance and extended infection. Antia et al. [3] have considered variant-transcending immunity as 
a basis for competition between variants, which can promote oscillatory behaviour, but this failed 
to induce sequential expression. Many other mathematical models of antigenic variation have been 
proposed and studied in the literature, but the discussion of their individuals merits and limitations 
is beyond the scope of this work. 

The model considered in this paper is a modification of the model proposed by Recker et al. |44] 
(to be referred to as Recker model), which postulates that in addition to a highly variant-specific 
immune response, the dynamics of each variant is also affected by cross-reactive immune responses 
against a set of epitopes not unique to this variant. This assumption implies that each antigenic 
variant experiences two types of immune responses: a long-lasting immune response against epitopes 
unique to it, and a transient immune response against epitopes that it shares with other variants. 
The main impact of this model lies in its ability to explain a sequential appearance of antigenic 
variants purely on the basis of cross-reactive inhibitory immune responses between variants sharing 
some of their epitopes, without the need to resort to variable switch rates or growth rates (see 
Gupta [27] for a discussion of several clinical studies in Ghana, Kenya and India, which support 
this theory). 

In the case of non-decaying long-lasting immune response, numerical simulations in the original 
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paper [33J showed that eventually all antigenic variants will be cleared by the immune system, with 
specific immune responses reaching protective levels preventing each of the variants from showing up 
again. Blyuss and Gupta [9] have demonstrated that the sequential appearance of parasitemia peaks 
during such immune clearance can be explained by the existence of a hypersurface of equilibria in 
the phase space of the system, with individual trajectories approaching this hypersurface and then 
being pushed away along stable/unstable manifolds of the saddle-centres lying on the hypersurface. 
They also numerically analysed robustness of synchronization between individual variants. Under 
assumption of perfect synchrony, when all variants are identical to each other, Recker and Gupta 
|46| have analysed peak dynamics and threshold of chronicity, while Mitchell and Carr [38J have 
considered the case of slowly decaying specific immune response. De Leenheer and Pilyugin |17] 
have replaced linear growth of antigenic variants in the original model by the logistic growth, 
and have studied the effects of various types of cross-reactivity on the dynamics, ranging from no 
cross-reactivity to partial and complete cross-immunty. 

It is known that time delay in the immune response can have a profound effect on the dynamics 
of host-parasite interactions and the host ability to eliminate infection. Several models have studied 
mathematically the effects of time delay on the immune dynamics and possible onset of oscillatory 
behaviour [131 l35| l36l I37j . For the Recker model, Mitchell and Carr have investigated the effect 
of time delay in the development of immune response in the case of complete synchrony between 
antigenic variants [38], and they have also investigated the appearance of synchronous and asyn- 
chronous oscillations [39J in the case of global coupling between variants (referred to as "perfect 
cross immunity" in |17j). 

In this paper, we use methods of equivariant bifurcation theory for delay differential equations 
to study the dynamics of a fully symmetric state in the Recker model. Stability analysis of the 
appropriate characteristic equation will show that under certain conditions on parameters, this 
state can undergo Hopf bifurcation, giving rise to different types of stable periodic solutions. The 
outline of this paper is as follows. In the next section we introduce the mathematical model 
of antigenic variation with time delay in the immune response and discuss its main properties. 
Section 3 discussed different steady states and derives the transcendental characteristic equation, 
which determines the stability of the fully symmetric state. In Sections 4 and 5 we analyse the 
case of small and arbitrary time delay, respectively, and find the boundary of Hopf bifurcation in 
terms of system parameters and the immune delay. Section 6 contains numerical simulations of the 
model and analysis of the symmetry properties of different periodic solutions. The paper concludes 
in Section 7 with the discussion of results and an outlook. 

2 Mathematical model 

In this section, the model of the immune response to malaria is presented together with some facts 
about the dynamics of this system. Following Recker et al. [33], we assume that each antigenic 
variant % consists of a single unique major epitope, that elicits a long-lived (specific) immune 
response, and also of several minor epitopes that are not unique to the variant. Assuming that all 
variants have the same net growth rate <f>, their temporal dynamics is described by the equation 



where a and a' denote the rates of variant destruction by the long-lasting immune response Zi and 
by the transient immune response W{, respectively, and index i spans all possible variants. 

It is known that the discovery of infection by immune receptors does not instantaneously lead to 
the development of the corresponding immune response [35]. To include this feature explicitly in the 
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model, we introduce time delay r as the time delay that elapses between changes in parasitemia and 
production of the corresponding immune effectors [381 [55] . For simplicity, it will be assumed that 
this time delay is the same for both specific and cross-reactive immune responses. The dynamics 
of the variant-specific immune response can be written in its simplest form as 



dz- 

= PVi(t ~ r) - fiZi(t), 



(2) 



where f3 is the proliferation rate, /i is the decay rate of the immune response, and r is the above- 
mentioned time delay in the immune response. Finally, the transient (cross-reactive) immune 
response can be described by a minor modification of the above equation ([2]): 



dwi 
~dT 



H Wi, 



(3) 



where the sum is taken over all variants sharing the epitopes with the variant j/j. We shall use 
the terms long-lasting and specific immune response interchangeably, likewise for transient and 
cross-reactive. 

The above system can be formalized with the help of an adjacency matrix T, whose entries 
Tij are equal to one if the variants i and j share some of their minor epitopes and equal to zero 
otherwise. Obviously, the matrix T is always a symmetric matrix. Prior to constructing this 
matrix it is important to introduce a certain ordering of the variants according to their epitopes. 
To illustrate this, suppose we have a system of two minor epitopes with two variants in the each 
epitope. In this case, the total number of variants is four, and they are enumerated as follows 



1 

2 
3 
4 



11 
12 
22 
21 



(4) 



The diagramme of interactions between these antigenic variants is shown in Fig. [TJ It is clear that 
for a system of m minor epitopes with rij variants in each epitope, the total number of variants is 
given by 



N 



n 



71; 



(5) 



After the ordering of variants has been fixed, it is straightforward to construct the connectivity 
matrix T of variant interactions. For the particular system of variants @, this matrix has the form 
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(6) 



For the rest of the paper we will concentrate on the case of two minor epitopes, but the results 
can be generalized to larger systems of antigenic variants. Using the connectivity matrix one can 
rewrite the system ([!])- ([3]) in a vector form 




a'w) 



-F(y,z,w) 



y{4>l N ~ az 
/3y{t - t) - iiz, 
(3'Ty(t - r) - /i'w 



(7) 
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Figure 1: Interaction of malaria variants in the case of two minor epitopes with two variants in 
each epitope. 

where y = (yi,y2, ■■■,Un) etc., Ijv denotes a vector of the length N with all components equal to 
one, and in the right-hand side of the first equation multiplication is taken to be entry-wise so 
that the output is a vector again. The above system has to be augmented by appropriate initial 
conditions, which are taken to be 

y(6) = i/>(9)>0,9e[-T,0], 
z(0) > 0,w > 0. 

with the history function ip(9) G C([— r, 0], R ), where C(\—t,0],M. n ) denotes the Banach space of 
continuous mappings from [— r, 0] into M. N equipped with the supremum norm \\tp\\ = sup_ T<e<0 \ip{9)\ 
for ifi G C([— r, 0], IR^), where | • | is the usual Euclidean norm on R . Using the same argument as 
in [9], it is possible to show that with these initial conditions, the system ([7|) is well-posed, i.e. its 
solutions remain non-negative for all t > 0. 

We will assume that cross-reactive immune responses develop at a slower rate than specific 
immune responses, have a shorter life time, and are less efficient in destroying the infection. This 
implies the following relations between the system parameters 

a < a, pi < fjf j < /3. (8) 

In terms of symmetry in the network of interactions between different antigenic variants, in the 
case of two minor epitopes with m variants in the first epitope and n variants in the second, the 
system ([7]) is equivariant with respect to the following symmetry group [3 [9] 

-p J S m x S n , vn 7^ n, , , 

| S m xS m xZ 2 , ra = n. 

Here, S m denotes the symmetric group of all permutations in a network of n nodes with an all-to-all 
coupling, and Z2 is the cyclic group of order 2, which corresponds to rotations by it. 

The above construction can be generalized in a straightforward way to a larger number of minor 
epitopes. System ([7]) provides an interesting example of a linear coupling, which for N > 4 does not 
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reduce to known symmetric configurations, such as diffusive, star or ali-to-all coupling |42j. A really 
important aspect is that two systems of antigenic variants with the same total number of variants 
iV may have different symmetry properties as described by the symmetry group T depending on m 
and n, such that N = ran. 



3 Steady states 

In the particular case of non-decaying specific immune response (/x = 0) and instantaneous immune 
response (r = 0), steady states of the system (JJJ) are not isolated but rather form an iV-dimensional 
hypersurface Hq = {(y,z, w) E M. N : y = w = Otv} in the phase space [9J. Linearization near each 
equilibrium on this hypersurface has the eigenvalues (—a') of multiplicity N, zero of multiplicity 
N, and the the rest of the spectrum is given by 

(f) — az\, <f) — az2, ...,</)— azN. 

This suggests that the hypersurface consists of saddles and stable nodes, and besides the original 
symmetry of the system it possesses an additional translational symmetry along the z axes. Fur- 
thermore, stability of equilibria on the hypersurface does not depend on the time delay r. The 
existence of this hypersurface of equilibria in the phase space leads to a particular behaviour of 
phase trajectories, which mimics the occurrence of sequential parasitimea peaks in the immune 
dynamics of malaria HI] . 

When fi > 0, the structure of the phase space of the system (J7|) and its steady states is drastically 
different. Now, the only symmetry present is the original symmetry T, and the hypersurface of 
equilibria Hq disintegrates into just two distinct points: the origin O, which is always a saddle, and 
the fully symmetric equilibrium 

E = (Y1 N ,Z1 N ,W1 N ), where 



a fifi' + a'n c fi'fi aj3fi' + a'n c f5'/j, afi/J + a'n c P'fj, 

Here n c is the total number of connections for each variant, which in the case of two minor epitopes 
with m variants in the first epitope and n variants in the second, is equal to n c = m + n — 1. It 
has been previously shown that in the absence of time delay, the fully symmetric equilibrium E 
can undergo Hopf bifurcation as the system parameters are varied 019]. It is worth noting that 
if one assumes all variants to be exactly the same, the system collapses onto a system with just 3 
dimensions, but in this case it is possible to show that the fully symmetric equilibrium is always 
stable for r = [46J and can have a Hopf bifurcation for t > [38 . Besides the origin and the 
fully symmetric steady state, system ([7]) has 2^ — 1 other steady states, all of which are unstable 

mm- 

In order to understand the structure of the solution that arises from the Hopf bifurcation of the 
fully symmetric steady state E, we concentrate on a specific connectivity matrix T corresponding 
to a particular example of a system of two epitopes with two variants in each epitope, as shown in 
Fig. [TJ In this case, the system (|7|) is equivariant under the action of a dihedral group D4, which is 
a symmetry group of a square [7J . We can use the subspaces associated with four one-dimensional 
irreducible representations of this group to perform an isotypic decomposition of the full phase 
space M 12 as follows [3 HI [53]: 



>12 



{(1, 1,1, !),(!, -1,1,-1), (1,0, -1,0), (0, 1,0, -l)} 3 . (11) 
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The characteristic matrix of the linearization of system ([7]) near a fully symmetric steady state E 
has the block form 



J(A,r) 




4 




where O4 and I4 are 4x4 zero and unit matrices, and T is the connectivity matrix ([6]). Rather 
than compute stability eigenvalues directly from this 12 x 12 matrix, we use isotypic decomposition 
(jllj) to rewrite this characteristic matrix in the block-diagonal form [53] 



A(A,r) 



/ M + 2N 


3 


O3 


O3 \ 


O3 


M -2N 


O3 


O3 


3 


3 


M 


O3 


V o 3 


O3 


O3 


M ) 



(12) 



M 



(13) 



where 

—a —a' \ 

-(X + fj) , N = 

-(a + mO/ 

Here, matrix M is associated with self-coupling, and N is associated with nearest-neighbour cou- 
pling. From the perspective of stability analysis, eigenvalues of the characteristic matrix A(A, r) 
are determined as the the roots of the corresponding characteristic equation 





where 



det A(A,t) = (det[A!(A,r)]) 2 • det[A 2 (A,r)] • det[A 3 (A,r)] = 0, 
Ai,2, 3 (A,r) = 



(14) 



-A 

(3e- Xr 



-a 



B 'l, 



-At 



2,3 C 



—a 

-(A + /i) 

-(a + mO 



with B[ 23 = (3', 3(3', -(3' for matrices M, M + 2N and M - 2N, respectively. Stability is now 
determined by the roots of the corresponding characteristic equation 



A 3 + (/x + //)A 2 + w'\ + \( a f3 + a'B')e' Xr + (a/3/i' + a ' B' 'fi)e~ XT 
This equation can be rewritten in the form 

A 3 + ,4A 2 + BX + CXe- Xr + De~ XT = 0, 

where 



0. 



(15) 
(16) 



A = n + y! , B = /z/x', C = aj3 + a'B', D = a/V + a'B'fi. 
When the immune response is instantaneous (r = 0), the above equation simplifies to 

A 3 + AX 2 + (B + C)X + D = 0, 

and one can use the Routh-Hurwitz criterion to deduce that in the cases B' = f3' and B' = 3/3', the 
fully symmetric steady state E never loses stability, as the eigenvalues remain in the left complex 
half plane for all values of the parameters. When B' = —f3', the steady state E can undergo Hopf 
bifurcation at 

, aPu + nil' (11 + n') 
a " = J^' ' (17) 
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thus giving rise to an odd periodic orbit, where variants 1 and 3 are synchronized and half a period 
out-of-phase with variants 2 and 4, i.e. each variant is tt out of phase with its nearest neighbours 
[TJ. Another possibility is for the steady state E to undergo a steady state bifurcation at 

, afifj,' 

provided fjf < afi\ij (at.fi — // 2 ), but due to restrictions on parameters flSJ), this cannot happen. 

Cubic quasi-polynomial equations similar to (fl~6j) have been previously studied in several applied 
contexts, such as, models of business cycles [H], testosterone secretion [HI HE], or neural networks 
with bidirectional associative memory [3§]. In each of those cases, analyses of the appropriate 
characteristic equation allowed one to find restrictions on system parameters and the time delay, 
which lead to the occurrence of Hopf bifurcation. 

We will use the results of equivariant bifurcation theory for delay differential equations to 
analyse symmetry properties of possible solutions arising at the Hopf bifurcation. While the effects 
of symmetry on the dynamics of Hopf bifurcation in systems without time delay have been known 
for quite a long time [Bl EH Ell EH], it is only in the last 10-15 years that these results have been 
adapted to the analysis of delay and functional differential equations. In a series of papers, Wu 
and co-workers extended the theory of equivariant Hopf bifurcation to systems with time delays 
and employed equivariant degree theory to study existence, multiplicity and global continuation of 
symmetric periodic solutions [29] [30], EU [55] . The results of this analysis have been subsequently 
applied to the studies of Hopf bifurcation in a number of symmetric models of coupled oscillators 
with delayed coupling pU US QH E6j [56] . 

The strategy now is to consider when the eigenvalues of the characteristic equation (|16p cross 
the imaginary axis with non-zero speed, giving rise to a Hopf bifurcation. One has to separately 
consider the cases B' = 3/3', B' = —fi', and B' = fi' , as these correspond to a Hopf bifurcation in 
the even, odd and V4 subspaces, respectively [53J. In the case of bifurcation in the even subspace, 
the periodic solution that appearing will be an even periodic orbit, which has the full original 
symmetry D4 and is characterized by all variants oscillating in perfect synchrony with each other. 
An odd solution, also called anti-phase solution, corresponds to a bifurcation in the odd subspace, 
and has variants 1 and 3 oscillating synchronously and half of a periodic out-of-phase with variants 
2 and 4. Finally, when the bifurcation takes place in the V4 subspace, this is known as the Hopf 
bifurcation with symmetry [53]; as is clear from (|14p . in this case, two pairs of complex conjugate 
eigenvalues simultaneously cross the imaginary axis, and this can give rise to a number of different 
periodic behaviours, including edge and vertex oscillations, as well as discrete travelling waves 
[53j . By finding the minimum value of time delay, at which Hopf bifurcation occurs for one of the 
three possible values of B' , one can identify the type of periodic solution that will appear at the 
corresponding Hopf bifurcation. 

4 Small immune delay 

We begin our analysis of stability by considering the case when the time delay in the development 
of immune response is small < r <C 1. In this case, one can write e~ Ar ~ 1 — At, which transforms 
the characteristic equation (]16p into a regular cubic equation 

A 3 + ai(r)A 2 + a 2 (r)A + a 3 = 0, (18) 
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where 



ai (r) = fi + fi' - r(aj3 + a'B'), 



a 2 (r) = /i/i' + a/3 + a'B' - r(a/V + a'B'fj,), (19) 
03 = a/3// + a'B' 11. 

While finding the roots of the cubic equation f|18jl is still too cumbersome, one can resort to the 
Routh-Hurwitz criterion to establish the conditions for stability of fully symmetric steady state E. 
These conditions are as follows 

a,i > 0, i = 1, 2, 3, 
(1\(12 — 03 > 0. 

For r small but different from zero, we can find that the Hopf bifurcation will occur when 



a\(r) > 0, 02(7") > 0, £13 > 0, and 
oi(t)o 2 (t) = a 3 , 



(20) 



which gives the critical value of the time delay as 



p- ypi-iQR 

= 2Q ' (21) 

with 

Q = (a/3// + a'B'/x)(a/3 + a'B'), 

P=(a(3 + a'B')(a/3 + a'B' + /x//) + (/i + //)(a/V + a'B'/i), 

R = a p [j, + a'B' y + ////(/i + //)• 
If we define a characteristic polynomial as 

5 (A, r) = A 3 + ai(r)A 2 + a 2 (r)A + a 3 , 

the characteristic equation (fT8|) at t = th becomes 

g(X,r H ) = A 3 + ai(r^)A 2 + a 2 (r^)A + ai(T H )a 2 (T H ) = 0. 

The eigenvalues of (I18p at r = t# can be readily found as 

Ai(rff) = -ai(r ff ) < 0, 

and 

where lv is the Hopf frequency. To establish the occurrence of a Hopf bifurcation, we need to show 
that Re(dX / dr)\ T=T[1 > 0. Since dg/dr = (dg/dr) + (dg/dX)(dX/dr) = 0, we have 

d\ _ dg idg _ -(a/3 + a'B')A 2 - oi(r)o 2 (r)A 
d7 ~ ~drl dX ~ 3A 2 + 2ai(r)A + a 2 (r) ' 

Evaluating this at r = Th gives 
dA(r 



fir 



[a 2 (Tff)(a/3 + a'B') - ai(T# )a 2 (Ttf) 3 / 2 z] [-a 2 (r#) - 01(7^)^/02(^)1] 
2a2(rff)[ai(rff) + a 2 (T H )} 
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Figure 2: Critical time delay th (|2ip at the Hopf bifurcation as a function of a'. Parameter values 
are P = 1, a = 1.2, /i = 0.1, = 0.8, (a) // = 0.1, (b) y! = 0.127, (c) y! = 0.18, (d) = 1.2. The 
colour corresponds to B' = j3' (blue), B' = 3/3' (green), and B' = —(3' (red). 

The real part can be found as 

"dRe(A(r)) 

When 5' = /3' or 5' = 3/3', this is true since all parameters are positive; for B' = —/3', this inequal- 
ity holds due to a' < a and j3' < /3, as required by ([8]). This implies that in all three cases, the 
eigenvalues A2,3 cross the imaginary axis at r = th with a positive velocity, which implies the ex- 
istence of a Hopf bifurcation at r = th- These finding can be summarized in the following theorem. 

Theorem 1. For sufficiently small values of the time delay 0<t<1, the fully symmetric steady 
state E is stable provided r < th as defined in h21\) and unstable for r > th . At r = th, this 
steady state undergoes Hopf bifurcation. If the minimal value of th corresponds to B' = 3/3', the 
bifurcating solution will be an even periodic orbit; if it corresponds to B' = —f3', the bifurcating 
solution will be an odd periodic orbit; if it corresponds to B' = (3 1 , the bifurcating periodic orbit will 
lie in the V4 subspace. If a' H as defined in |i7p satisfies a' H < a, then for a' = a' H , the steady state 
E undergoes Hopf bifurcation to an odd periodic orbit for zero time delay r. 



Q2(th)[o? (t h ) + a/3 + a'B'} 
2[ai(rff) + a 2 (T H )] 



> 0. 



(23) 
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Figure [2] illustrates how the critical time delay th at the Hopf bifurcation varies with the rate a' 
of variant destruction by transient immune response, and also with the decay rate of the transient 
immune response //. The type of bifurcating periodic orbit is determined by which of the matrices 
M, M ± N will become unstable first as the time delay r increases. One can observe that when 
the death rate of cross-reactive immune response \j! is sufficiently close to that of the specific 
immune response fx, i.e. when the cross-reactive immune response is sufficiently long-lasting, the 
only possible periodic solution the steady state E can bifurcate to is the even solution (see Fig. (a)), 
in which all antigenic variants are behaving identically to each other. As the lifetime of transient 
immune response gets shorter (i.e. the rate \J increases), there appears a range of a 1 values shown 
in Figs, (b) and (c), for which the steady state E can also bifurcate to an odd periodic orbit, in 
which variants 1 and 3 are synchronized and half-a-period out of phase with variants 2 and 4. As 
\J increases further, the possibility of bifurcating into an even solution completely disappears (see 
Fig. (d)), and the only remaining possibility is a bifurcation to an odd solution for all admissible 
values of a'. It is worth noting that for higher values of //, the range of possible values of a' for 
which Hopf bifurcation can occur, is bounded by a' H given in (I17j) . for which the steady state E 
bifurcates to an odd periodic orbit for r = 0. 

5 The case of general delay 

In the situation when the immune delay r is not small, the approximation used in the previous 
section is not valid, and one has to consider the full equation (TIB]) . Before proceeding to the 
analysis of the full characteristic equation, it is worth mentioning a general result of Hale [28], which 
guarantees absolute stability of a delay system - this is a case when the real parts of all eigenvalues 
remain negative for all values of time delay, i.e. effectively independent of the delay. This is achieved 
when the corresponding ODE system is asymptotically stable, and the characteristic equation has 
no purely imaginary roots. 

Now we consider the characteristic equation (|16|) and analyse its roots in order to identify 
possible parameter regimes when the fully symmetric steady state E can lose its stability. First of 
all, due to biological restrictions on system parameters, the coefficient D in the equation (|16p is 
always positive, and this implies that A = is a not a root of the equation. Therefore, the only way 
how the steady state E can lose its stability is through a Hopf bifurcation at some r = tq, when a 
pair of eigenvalues of the equation (|16p crosses the imaginary axis. Let us assume that when r = 0, 
the steady state E is stable, which is always the case for B' = (3' and B' = 3/3', and is also true 
for B' = —fi 1 provided a' < a' H , where a' H is given in (fTTj) . In order to find out when the Hopf 
bifurcation can occur as the time delay r increases from zero, we look for solutions of equation (|16p 
in the form A = ioo (u > 0). Such a solution would be a root of (|16|) if and only if to satisfies 



Squaring and adding these two equations yields a single polynomial equation for the Hopf frequency 
to: 




Separating the real and imaginary parts, we have 




(24) 




(25) 



Let x = to 2 , and also 



ci = A 2 - 2B, c 2 = B 2 - C 2 , c 3 = —D 2 . 



(26) 
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Now, introducing function f(x) as 



f(x) = x 3 + c\x 2 + c 2 x + c 3 , 



(27) 



we can rewrite equation (|25p as 



f(x) = X 3 + C\X 2 + C2X + C3 = 0. 



(28) 



Since /(0) = 03 = — D 2 < 0, and lim x _ i . +00 /(x) = +00, we conclude that equation h{z) always has 
at least one real positive root. Without loss of generality, suppose equation ([28]) has three distinct 
positive roots, denoted y xi, X2 and X3, respectively. Then, correspondingly, equation ([25]) also has 
three positive roots 

= y/xi, UJ 2 = y/x2, W3 = y/x3. 

From (1241) one can find 



COS OJ^T 



u 2 (AD -CB + Cuf, 
D 2 + C 2 co 2 



, k = 1,2,3. 



If we denote 



» 



arccos 



(AD - gff + go; 2 ) 
L> 2 + C 2 w 2 



+ 27rn 



1,2,3, n = 0,l,2, 



(29) 



00, 



(30) 



then ±iuik is a pair of purely imaginary roots of ([To]) with r = . Since lim n _ s . 0O t^" - -* 
A; = 1, 2, 3, we can define 

(0) ■ r (0)-, 

Using the results of [48J, we can conclude that all roots of the characteristic equation (|16p have 
negative real parts when r £ [0, To). By construction, it follows that /'(wq) > 0, and this implies 
that at t = tq, one has ±iuo as a pair of simple purely imaginary roots of equation (fT6|) . see [37] 
for details of the proof. The next step is that show that, in fact, 



dRe(A(r)) 
dr 



> 0. 



To do this, we differentiate both sides of equation (I16p with respect to r, which yields 

dX{r) {CX 2 + DX) e~ Xr 

dr ~ 3A 2 + 2AX + B + {C - CtX- TD)e~ Xr ' 

Evaluating the real part of this expression at r = tq and substituting A = gives 



rfRe(A(r)) 
d7 



lu; 



[3cu$ + 2lu 2 {A 2 - 2B) + B 2 -C 2 



A 



T=TQ 



where 



(31) 



A = (— 3wq + B — tqAloq + C cos wqtq) + (2Ajq + tqo;o(-B — Wq) — C sinwo^o)' 
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(d) 




Figure 3: Critical time delay tq fl30f) at the Hopf bifurcation as a function of a' . Parameter values 
are /3 = 0.2, a = 0.3, £t = 0.1, ft = 0.8, (a) y! = 0.1, (b) y! = 0.24, (c) // = 0.255, (f) // = 0.7. 
The colour corresponds to B' = (3' (blue), B' = 3(3' (green), and B' = —j3' (red). 

Using the definition of function / from ([27]) with the coefficients given in (I26p , we can alternatively 
write 



dRe(A(r)) 



dr 



A 



(32) 



Since /(0) < and ujq is the smallest positive root of (j25|) . it follows that 

unless ujq is a double root, in which case we take u)q as the next root. Therefore, one concludes that 



dRe(A(r)) 



A 



>0, 



which, in the light of Hopf theorem, implies that at r = To, the fully symmetric steady state un- 
dergoes a Hopf bifurcation. We summarize these results as follows. 



Theorem 2. The fully symmetric steady state E is stable for r < tq, where tq is defined in i30\) . 
unstable for t > tq, and undergoes a Hopf bifurcation at r = tq. If the minimal value of To cor- 
responds to B' = 3/3' , the bifurcating solution will be an even periodic orbit; if it corresponds to 
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B' = —ft' , the bifurcating solution will be an odd periodic orbit; if it corresponds to B' = f3' , the 
bifurcating periodic orbit will lie in the V4 subspace. 

Figure [3] illustrates how the time delay affects possible types of bifurcating solutions. Similar to 
the case of small time delay, the steady state E bifurcates primarily into an even periodic orbit for 
smaller values of \j! and into an odd periodic orbit for higher values of /j,'. The main difference from 
the case of small time delay is that now the range of admissible values a' is bounded by a rather 
than the value of a' H at the Hopf bifurcation for r = 0. There are two main conclusions to be drawn 
from this Figure. The first one is that the computations of tq suggest that the state E can only 
bifurcate to either even or odd periodic orbit, and the bifurcation into a subspace V4 corresponding 
to two pairs of complex conjugate eigenvalues crossing the imaginary axis simultaneously does not 
happen, thus excluding the occurrence of vertex or edge oscillations, as well as discrete travelling 
waves [S3]. This, however, does not preclude them from appearing in the dynamics altogether, as 
they can arise as solutions bifurcating from the even/odd periodic orbits, as will be shown in the 
next section. The second conclusion which follows from Fig. [3] is that as the efficiency of cross- 
reactive immune response a' increases, this leads to a decrease in To, thus leading to the onset of 
sustained oscillations for faster immune responses. 

6 Numerical simulations 

In the previous sections we established that the fully symmetric steady state E can undergo Hopf 
bifurcation, giving rise to two different types of solutions: a symmetric solution, in which all 
variants are oscillating identically, and an odd periodic orbit having variants 1 and 3 oscillating 
synchronously and half a period out-of-phase with variants 2 and 4. To understand evolution of 
these solutions as the time delay increases, we have performed a number of numerical simulations of 
system ([7|) , which are shown in Fig. HI One can observe that for sufficiently small immune delay, the 
fully symmetric steady state E is stable, see Fig. (a). As the time delay r exceeds its critical value 
(I30|) at the Hopf boundary, this state becomes unstable and gives rise to a stable fully symmetric 
periodic orbit, as demonstrated in Fig. (b). As the time delay increases, this solution acquires 
subharmonics, as shown in Fig. (c), and it eventually becomes chaotic, see Fig. (d). For other 
values of parameters, the fully symmetric steady state E bifurcates into an odd periodic solution 
(see Fig. (e)), which for higher values of r transforms into a discrete travelling wave, where each 
of the variants is quarter of period out-of-phase with its neighbours on the diagramme (p}. In fact, 
as the time delay is varied, it is possible to observe other stable periodic solutions with different 
phase shifts between antigenic variants. 

The important difference from the model without time delay is that the only possibility for 
r = is a bifurcation into an odd periodic orbit [7], while for r > there is also a possibility of the 
steady state E bifurcating to a stable fully symmetric periodic orbit. In this case, Hopf bifurcation 
does not lead to the breakdown of the original D4 symmetry. Also, we note an important difference 
from the globally coupled system with delayed immune response. When considering a network 
of antigenic variants with an all-to-all coupling, Carr and Mitchell [39] have shown that for the 
majority of parameter values, anti-phase Hopf bifurcation eventually leads to the behaviour that 
appears chaotic, while the simulations shown in Fig. 2] suggest that when not all antigenic variants 
are related to each other, the system is able to support a number of different out-of-phase solutions 
without going into chaotic regime. 

To better understand the structure of different types of periodic behaviour in the model, one can 
use the so-called H/K Theorem, which takes into account individual spatial and spatio-temporal 
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Figure 4: Numerical solution of the system ([7j) with a connectivity matrix ([6]). Parameter values 
are <f> = 1, (3 = 0.2, a = 0.3, \i = 0.1, /?' = 0.8. Different colours correspond to different antigenic 
variants, (a) Stable steady state E (// = 0.26, a' = 0.1, r = 1). (b) Even periodic orbit (// = 0.26, 
a' = 0.1, t = 1.26). (c) Quasi-periodic even orbit (// = 0.26, a' = 0.1, r = 1.27). (d) Chaotic 
solution (// = 0.26, a' = 0.1, r = 1.5). (e) Odd periodic orbit {jl = 0.4, a' = 0.2, r = 1.2). (f) 
Discrete travelling wave (// = 0.7, oi = 0.2, r = 1.8). 

symmetries of the solutions [121 125j . To use this method, we note that due to T-equivariance of the 
system ([7]) and uniqueness of its solutions, it follows that for any T-periodic solutions x(t) and any 
element 7 G T of the group, one can write 

7x(i) = x(t - 9), 

for some phase shift 9 G S 1 = R/Z = [0,T). The pair (7, 0) is called a spatio-temporal symmetry 
of the solution x(t), and the collection of all spatio-temporal symmetries of x(t) forms a subgroup 
A C f x S 1 . One can identify A with a pair of subgroups, H and K, such that K C H C T. We 
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also define 



H = {7 G r : ^y{x(t)} = {x(t)}} spatio-temporal symmetries, 
K = {7 G T : 7x(t) = x(t) Vt} spatial symmetries. 

Here, K consists of the symmetries that fix x(t) at each point in time, while H consists of the 
symmetries that fix the entire trajectory. Under some generic assumptions on H and K, the H/K 
Theorem states that periodic states have spatio-temporal symmetry group pairs (H, K) only if H/K 
is cyclic, and K is an isotropy subgroup |12l 125] . The H/K Theorem was originally derived in the 
context of equivariant dynamical systems by Buono and Golubitsky [12J, and it has subsequently 
been used to classify various types of periodic behaviours in systems with symmetry that arise in a 
number of contexts, from speciation [50J to animal gaits [43J and vestibular system of vertebrates 
[22]. 

In the case of D4 symmetry group acting on four elements, there are eleven pairs of subgroups H 
and K satisfying the above requirements [25J. While periodic solutions corresponding to each such 
pair can exist in a general setup, the above theorem does not guarantee their existence or stability 
in a particular system, see discussion in |25| . Therefore, we use numerical simulations shown in in 
Fig- HI to identify specific types of periodic solutions and their spatio-temporal symmetries that can 
be found in the system (J7|). Solutions shown in Figs, (b) and (c) have the full original symmetry of 
the system, and therefore are characterized by a pair (H, K) = (D4, D4). The solution shown in Fig. 
(d) is chaotic and does not have any of the symmetries of the original system. Figure (e) illustrates 
an odd periodic orbit with a symmetry (H, K) = (D4, DlJ), where K = Dr? is an isotropy subgroup 
associated with reflections along the diagonals of the square. Finally, the solution shown in Fig. (f) 
is a discrete travelling wave with the symmetry (H, K) = (Z4, 1), also known as a "splay state" [52], 
"periodic travelling (or rotating) wave" [6], or "ponies on a merry-go-round" or POMs [5] in the 
studies of systems of coupled oscillators. In this dynamical regime all variants appear sequentially 
one after another along the diagramme in Fig. ([I]) with quarter of a period difference between 
two neighbouring variants. From the perspective of equvariant bifurcation theory, this solution is 
generic since the group Z n is always one of the subgroups of the D n group for the ring coupling, 
or the S n group for an all-to-all coupling, and its existence has already been extensively studied 
[5| 1231 [2^]. From the immunological point of view, this is an extremely important observation that 
effectively such solution, which immunologically represents sequential appearance of parasitemia 
peaks corresponding to different antigenic variants, owes its existence not to the individual dynamics 
of antigenic variants, but rather to the particular symmetric nature of cross-reactive interactions 
between them. This immunological genericity ensures that the same conclusions hold for a wide 
variety of immune interactions between human host and parasites, which use antigenic variation 
as a mechanism of immune escape, as illustrated, for instance, by malaria, African Trypanosomes, 
several members of Neisseria family (N. meningitidis and iV. gonorrhoeae), Borrelia hermsii etc. 

[23 [53]. 

7 Discussion 

In this paper we have used methods of equivariant bifurcation theory to understand the effects 
of immune delay on the dynamics in a model of antigenic variation in malaria. In the simplest 
case of two epitopes with two variants in each epitope, the system is equivariant with respect to 
a D4 symmetry group of the square. Using isotypic decomposition of the phase space based on 
the irreducible representations of this symmetry group has allowed us to find critical value of the 
time delay at the boundary of Hopf bifurcation of the fully symmetric steady state in terms of 
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system parameters. We have identified even and odd periodic solutions that can arise at Hopf 
bifurcation, and also performed numerical simulations of the full system to illustrate other types 
of dynamical behaviours that can be exhibited by the model. These have been classified in terms 
of their spatial and spatio-temporal symmetries using the H/K Theorem. Our analysis suggests 
that as the efficiency of the cross-reactive immune response increases, the critical value of the 
time delay at the Hopf bifurcation decreases, which is similar to earlier studied cases of complete 
synchrony [46 [ [38] or global coupling between antigenic variants [39]. At the same time, unlike 
these , out system is able to support a range of stable phase shift solutions without developing 
chaotic dynamics. 

When applying the results of our analysis to the studies of realistic models of antigenic vari- 
ation, one of the important consideration that have to be taken into account is the fact that in 
reality systems of antigenic variants do not always fully preserve the symmetry assumed in the 
mathematical models. In the context of modelling immune interactions between distinct antigenic 
variants, this means that not all variants cross-react with each other in exactly the same quantita- 
tive manner. Despite this limitation, due to the normal hyperbolicity, which is a generic property 
in such models, main phenomena associated with the symmetric model survive under perturba- 
tions, including symmetry-breaking perturbations. The discussion of this issue in the context of 
modelling sympatric speciation using a symmetric model can be found in [25] . 




The results presented in this paper are quite generic, and the conclusions we obtained are valid 
for a wide range of mathematical models of antigenic variation. In fact, they are applicable to 
the analysis of within-host dynamics of any parasite, which exhibits similar qualitative features 
of immune interactions based on the degree of relatedness between its antigenic variants. The 
significance of this lies in the possibility to classify expected dynamical regimes of behaviour using 
very generic assumption regarding immune interactions, and they will still hold true, provided the 
actual system preserves the underlying symmetries. 

There are several ways in which the analysis in this paper can be further improved to achieve an 
even more realistic representation of the dynamics of antigenic variation. One of the assumptions 
in our analysis is that that the degree of cross-reactivity between antigenic variants does not vary 
with the number of epitopes they share. It is straightforward, however, to introduce antigenic 
distance between antigenic variants in a manner similar to the Hamming distance pQ [45] . Such 
a modification would not alter the topology of the network of immune interactions, but rather it 
would assign different weights to connections between different antigenic variants in such a net- 
work. Another modelling issue concerns the way how the time delay in the immune response can 
be represented mathematically in the most realistic manner. It would be beneficial to investigate 
the dynamics of antigenic variation under the influence of a more general distributed time delay, 
which is known to cause both destabilization of steady states, and also suppression of oscillations 
[H El [32~1 133j . This would provide better insights into how the efficiency of developing and main- 
taining immune response affects the within-host dynamics of parasites with antigenic variation. 
Alternatively, one can introduce different time delays for the development of specific and long- 
lasting immune responses and analyse how the difference in these timescales affects the overall 
stability and dynamics. Systematic analysis of the effects of the introduction of antigenic distances 
and different /distributed delays in the immune response on the dynamics of antigenic variation is 
the subject of further study. 
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